This tutorial covers how to run, inspect, and report a linear model in R. For the report portion, we will cover some key features of dynamic reporting in Quarto and how to write and render professionally formatted documents.
In this section, we will walk through the process of fitting, comparing, and reporting hierarchical linear models in R. This is not a statistics tutorial, so there will be minimal detail about how to understand or interpret the output of these commands. Instead, refer to Prof Andy Field’s discovr tutorials, which are the primary teaching materials for the same content in UG teaching. All of the code and interpretation in the following section is from discovr_08, the GLM.
There are two key goals for this linear model walkthrough:
discovr tutorialsOf the two, the first is more important. I’d strongly recommend you open the relevant discovr tutorial and skim through it as you go so you’re familiar with what it contains. However, the following sections of this tutorial will present the same code and information with very abbreviated text, to serve as a quick reference.
You can also have them both open at once and refer to each!
discovr tutorials
Prof Andy Field’s discovr tutorials provide detailed walkthroughs of both the R code and the statistical concepts of a variety of statistical analyses. They are a good place to look first to understand what your UG supervisees or advisees have been taught on a particular topic.
The tutorials are built in {learnr}, an interactive platform for learning and running R code. So, unlike the tutorial you’re currently reading, they must be run inside an R session. We have already installed all of the tutorials in your Posit Cloud workspace.
To open a tutorial, open any project and click on the Tutorial tab in the Environment pane. You can run any tutorial from here, but the relevant one for the linear modelling we are working on now is discovr_08, “the GLM”. Scroll down to this tutorial and click the “Start Tutorial ▶️” button to load the tutorial.
Because discovr tutorials run within R, you don’t need to use any external documents; you can write and run R code within the tutorial itself. However, I strongly recommend that whenever you work with these tutorials, you write and run your code in a separate document, otherwise you will have no record of the code and output.
Today’s data is a truncated version of the TeachingRatings dataset from the {AER} package. You can load the codebook in the Help viewer by running the following in the Console:
?AER::TeachingRatingsNow that we have some data, we can fit our first model with a single predictor. We will do this with the very hardworking lm() function in R, standing for “linear model”.
The lm() function has a lot of options (as you might expect, given the versatility and ubiquity of linear models!), but its basic format to fit a linear model with a single predictor is very simple:
lm(outcome ~ predictor, data = our_dataset)In this function, outcome ~ predictor is a formula expressing a (simplified version of) the linear model equation. Here, outcome is the name of the variable in our_dataset that contains the outcome or dependent variable y, and predictor is the name of the variable that contains the predictor or independent variable x.
Now we have a new object that contains all the information about our model. We could simply call the name of this object, but the output doesn’t tell us anything besides the actual value of the b estimates. Instead, we’ll use some useful functions from the {tidyverse} package {broom} to get the information we need.
To get some common measures of model fit, we can use the function broom::glance(). The output includes \(R^2\), adjusted \(R^2\), and F and accompanying statistics in comparison to the null model (the mean of the outcome alone).
Helpfully, broom::glance() (and many other {tidyverse} functions) output tibbles, which means we can work with them as we already know how to do. In Tutorials 4 and 5, we’ll also learn more about changing and filtering tibbles that will make this very more useful. For now, we can simply note the information we get out of this function for future use.
To get information about the b estimates for individual predictors, we can use the function broom::tidy(). We could run this function without any other information, as we did with glance() above, but we’ll change one argument here in order to get 95% confidence intervals for b in the output as well.
Next, we may want to test the addition of further predictors in the model. We can then compare the more complex multi-predictor model to the single-predictor model.
The models we’re describing here contain only additions, not interactions. So, you may be wondering, “What if I want to model more complex relationships between predictors?” If that’s what you’re trying to do now, you may want to jump ahead to discovr_10 for moderation and mediation.
Otherwise, we will get there in this course eventually!
Now we have two models, one simpler with only a single predictor, and the other with two predictors. We might next want to test whether the more complex, two-predictor model is a significant improvement over the simpler model, in order to decide which model to retain. We can do this with the anova() function1 to compare the two models.
The anova() function will only work for model comparison for particular models.
The models must be hierarchical. That is, the more complex model(s) must contain all of the predictor(s) present in the less complex model(s).
All models must be fit to the same dataset. If, for example, your first predictor has no missing values, but your second predictor had one, the model with using the second predictor would have been fit to a dataset of a different size than the model using only the first, and the anova() function will throw an error to this effect.
If you encounter this issue, you may need to inspect and clean your dataset before you proceed with analysis!
The output from many base-R {stats} functions, like anova(), include a line labeled Signif. codes that provide a key for understanding the notation given for significance levels in p-values.
Reading the key from left to right, we can see that a result is given three asterisks (***) when the p-value is between 0 and .001; two asterisks between .001 and .01; and so on.
This can be a useful visual check, especially because p-values that are very, very small are frequently expressed in scientific notation, which can make them more difficult to spot.
The F-test is significant, indicating that the more complex two-predictor model is a significant improvement over the one-predictor model, so we will proceed with the two-predictor model.
You may have noticed that the output we’ve seen so for only contains unstandardised model parameter estimates. If we want standardised Bs expressing the relationship between predictor(s) and outcome in standard deviation units, we can make use of the model_parameters() function from the {parameters} package to standardise our bs.
At Sussex, we teach a range of model diagnostics and robust model sensitivity tests in order to test model assumptions. We will look briefly at each of these in turn.
Remember, there are more examples and longer explanations in the discovr_08 tutorial!
To begin, we can generate nicely formatted, customisable residual plots using the function ggplot2::autoplot(). However, it is essential to load the {ggfortify} package in order for this to work!
Load the {ggfortify} package and use the autoplot() function to generate residual plots for eval_full_lm. Set the which argument to c(1, 3, 2, 4).
library(ggfortify)Warning: package 'ggfortify' was built under R version 4.3.1Loading required package: ggplot2Warning: package 'ggplot2' was built under R version 4.3.1ggplot2::autoplot(eval_full_lm, which = c(1, 3, 2, 4))
If you’re wondering what’s up with which, the plot() function that autoplot() is based on has a total of six plots available. Here I’ve chosen the two residual plots, the normal Q-Q and the Cook’s distance plot. Plots 5 and 6 have to do with leverage, which we don’t teach at UG level. To see them, simply add them to the which = argument.
If you want to customise the theme or look of these plots further, they are built with {ggplot2} so you can add or change anything about them using that package. We will come round to a detailed exploration of data visualisations with {ggplot2} in Tutorial 06.
In the UG core statistics modules at Sussex, we teach that normality is the least important of the assumptions of the linear model - the most important being additivity and linearity - so we do not generally worry too much about normally distributed standardised residuals, especially in large sample sizes. However, we do teach them how to evaluate the proportion of standardised residuals with values above $$1.96 (approximately 5%), $$2.56 (approximately 1%), and above $$3 (likely an outlier). For details on how to use broom::augment() to obtain standardised residuals and other model diagnostic measures like Cook’s distance, see the discovr_08 tutorial.
At UG level, we teach robust models for two purposes:
Our first robust model re-estimates the parameter estimates using robust techniques withe the robust::lmRob() function. We can then compare the robust parameter estimates to the unadjusted ones obtained with lm() to find out if our estimates were biased, and report the robust parameter estimates if they were.
To test and adjust for heteroscedastic residuals, we can re-estimate the standard error using a robust method. To do this, we’ll use the parameters::model_parameters() function with the argument vcov = "HC4" as recommended in the discovr tutorial. To do this, use the unadjusted model object eval_full_lm as the first argument.
If we had a small sample size, a final option would be to bootstrap the confidence intervals. To do this, we will again use parameters::model_parameters(), but this time with bootstrap = TRUE.
Note: Sample size is not an issue with this dataset (N = 463).
Quarto documents are a mix of text and code, the next generation of R Markdown documents. For our purposes, we will be using R within Quarto, but Quarto documents support the integration of many different coding languages, including Python, Julia, and Observable. If you’ve previously used Rmd (RMarkdown), Quarto is backwards-compatible and will able to render most documents with no issues.
The official Quarto Guide is extensive, detailed, and extremely helpful. It’s always the best first stop for any questions you have about using Quarto. Quarto also offers detailed tutorials.
This quick-reference to Markdown formatting is particularly helpful.
Prof Andy Field has also recorded a series of video guides to using Quarto that are used in UG teaching.
To get some hands-on practice working with Quarto, we will create a new Quarto document from scratch. We will use the analysis code we’ve already written to create a nicely formatted report.
By default, your new Quarto document will already have some settings and content to demonstrate how it works. Most importantly, you can see that there are three main types of information in this document:
---s, which contains information about how the document will be rendered
A new feature with Quarto is Visual mode, which is very like using Word or a similar word processing programme. All formatting can be applied and previewed in the document itself, using keyboard shortcuts or the familiar formatting buttons along the top toolbar.
Visual mode also has an “insert anything” shortcut, /, that allows you to quickly insert elements into your document. If you type / at the start of a new line (or Ctrl/Cmd + / otherwise), a drop-down list of possible elements will appear, which you can navigate by scrolling or search by typing.
Alternatively, you can toggle to Source mode by clicking the “Source” button at the very top left of the document pane. In Source mode you can edit the document using Markdown formatting. This allows more fine-tuned control of elements, but automatic formatting and the / shortcut don’t work in this mode.
Which you use is entirely personal preference - I toggle regularly between them depending on what I’m trying to do, so pick whichever works for you.
To prepare for our analysis, we will need a place to load packages and read in the data. Any executable code - that is, code that you want to run and do something - must be written in a code chunk and NOT in the body text2.
There are several ways to insert a new code chunk.
/ and then Enter (since “R Code Chunk” is the first option)```{r} ```To begin, we’ll use headings to map out our document. Properly formatted headings are strongly recommended for any documents you write, for a variety of reasons:
There are several ways to insert a new heading.
### creates a level 3 heading)Any new-line text will automatically be plain text. This can be anything you like, although in our current document, we are writing a mock-formal results section.
How to format body text depends on the mode you are in.
That last sentence should report which of the two models was better, based on the result of our F-test, rather than question marks. We could produce the output of this test, read it ourselves with our very own eyes/ears, and then type out the results by hand…but that is definitely not what we are going to do! Instead, we’ll look at a couple options for reporting the results dynamically. (We’ll come back to the model R2s in a moment.)
Our first option is to use inline code to report the numbers. In order to do this, we first need to know what information we have available in the test output, so we can make use of it.
This section will make extensive use of $ subsetting, which was covered in Tutorial 02, and [] subsetting, which was covered in Tutorial 01.
The tibble we get here is quite different from the “Analysis of Variance Table” text output that we saw earlier when we printed out the results of the anova() function on its own. We now have all the values we need to report these results in a conveniently subsettable tibble with nice R-friendly names.
So, let’s have a go getting at some of those values.
We now have a bit of R code that produces a single number that we want to report in the text3. The issue is that code chunks contain code, and body text contains text, but we would like the code that we’ve written to print out its value in the text!
The solution is inline code, a small bit of R code written in the body text (“inline”) that will, when the document is rendered, automatically insert the right number. Inline code is written as follows, in the text: `r some_r_code`
For our reporting, we will need to insert the inline code in exactly the spot where we would like the output of the code to appear. This particular bit of code produces the F-statistic, so we can replace the “???”s in our reporting where the value of the F-statistic should go:
An ANOVA comparing the models indicated a significant improvement in model fit for the second model compared to the first (F(???, ???) =
`r round(tidy_f$statistic[2], 2)`, p = ???).
You can check whether this has worked in two ways:
As helpful as inline code is (and I would recommend reporting all values dynamically/automatically wherever possible, so it is very useful!), you may have noticed that there was a lot of repetitive typing that also made the text itself quite difficult to read, as well lots of opportunities make typos or mistakes. Surely there’s a simpler way to do this sort of thing!
There are, in fact, many simpler ways to do common tasks like this, which take advantage of the fact that an object created by a particular function will always have the same structure. One option is to make further use of the {papaja} package, which is designed for just this purpose.
Next, we will jump to the “Final Model” heading and have a look at how to turn our final model output into a nicely formatted table. Once again, {papaja} provides a quick and beautiful solution for reporting, so let’s use it again.
As anyone who has had to create a long document with lots of tables and figures knows, keeping track of the numbering is a huge pain, especially when, for instance, a reviewer asks you to add something partway through and then everything has to be renumbered4.
The good news is that Quarto can take care of figure and table numbering automatically. There are two steps to this:
fig- for figures and tbl- for tables.@.After all this work to analyse, interpret, and report, it’s finally time to produce the final document. The process of turning a Quarto document into some output format - including running all the code and applying all the formatting - is called rendering (previously knitting with RMarkdown).
At the moment, your report may not be as clean as we’d like it to be: there are likely messages from R floating around and code all over the place, whereas for a formal report, we of course only want to show the final output. This is where the YAML header comes in.
The default YAML header contains only a few explicit settings, and if you haven’t changed anything, probably looks like this:
---
title: "Untitled"
format: html
editor: visual
---By adding options to the YAML header, we can determine how the document as a whole is rendered. Here are the common ones that I use on the regular:
---
title: "Linear Model Report"
format:
html:
1 toc: true
editor: visual
2self-contained: true
3execute:
4 echo: false
5 warning: false
6 message: false
---The requirements for each document will change depending on its purpose, so you will likely find that
The Quarto help documentation, as usual, has a complete list of YAML options, including how to set default behaviour for figures and other settings.
Global options apply to the entire document, but you may want to change these settings for individual code chunks to override the default settings in the YAML.
For example, you may have a code chunk containing some processing code that you used to view and clean your data. If your global execution option is set to echo: false, the code output would still appear, although the code itself would be hidden. If, for this particular code chunk, you don’t want the output to appear either, you can override the global option with a local option for that chunk only.
Local code chunk options appear as “hashpipe” (#|) comments, as we have seen earlier with labels. They use the same syntax as YAML, but the settings only apply to individual code chunks. Any settings that aren’t explicitly changed within the chunk are inherited from the YAML settings.
For this example, we could set a local code chunk option include: false which will prevent the output from appearing in the document.6
```{r}
#| include: false
some_code_doing_cleaning_and_processing
```For these tutorials, we will generally stick to HTML, as it’s the most painless of the rendering options. However, you will likely find yourself wanting to produce some other type of document, which you can easily7 do from the same Quarto document.
To render to a different format, change the YAML format: setting to a different output.
As per, the Quarto guide on output formats has all the information you need!
And there we have it! You now have a complete example of a linear model report, rendered into both HTML and Word, to refer to. It’s amazing how much you’re able to do after just a few weeks!
This is the end of the FundRmentals section of the course. If you’re so inclined, we’ll see you in the Essentials section, which will cover data wrangling and cleaning, and running and reporting many more statistical analyses.
Somewhat confusingly, this is not the function we teach UGs to perform ANOVAs! See discovr tutorials 11, 12, 15, and 16 for a detailed guide through the afex package for running linear models with categorical predictors.↩︎
Mostly true, except for inline code!↩︎
Generally, inline code should only ever produce a single value, otherwise the formatting can get interesting. This value, however, could longer than a single number, if you want to get creative with your dynamic reporting!↩︎
It seems obvious that this should be R instead of K, but remember this made perfect sense when it was called “knitting”!↩︎
For writing tutorials, I make extensive use of eval: false, which includes the code in the output but does not attempt to run the code. This allows me to write all kinds of nonsense code without R getting stroppy and throwing errors all of the place!↩︎
Depends crucially on what you consider to be “easy”, especially when dealing with PDFs!↩︎
